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We study stationary states in a diluted asymmetric (kinetic) Ising model. We apply the recently 
introduced dynamic cavity method to compute magnetizations of these stationary states. Depend- 
ing on the update rule, different versions of the dynamic cavity method apply. We here study 
synchronous updates and random sequential updates, and compare local properties computed by 
the dynamic cavity method to numerical simulations. Using both types of updates, the dynamic 
cavity method is highly accurate at high enough temperatures. At low enough temperatures, for 
sequential updates the dynamic cavity method tends to a fixed point, but which does not agree 
with numerical simulations, while for parallel updates, the dynamic cavity method may display 
oscillatory behavior. When it converges and is accurate, the dynamic cavity method offers a huge 
speed-up compared to Monte Carlo, particularly for large systems. 
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I. 'INTRODUCTION 



Stationary states of classical equilibrium systems are described by the Boltzmann-Gibbs measure. In complex 
systems, the exact computation of even local properties (marginals) is not feasible, and perturbative methods or 
other approximations therefore have to be used. Much attention has over the last decade been given to the Belief 
Propagation (BP) or Bethe-Peierls ansatz class of approximations of the marginals, which are exact if the underlying 
graph of interactions is a tree, and generally expected to be accurate if the underlying graph is locally tree- like [171 [33]. 

In contrast to equilibrium systems, non-equilibrium systems do not admit a similar, universal, description of their 
stationary states. We here take a small step in this direction. We show that the recently introduced dynamic cavity 
method, essentially a Bethe-Peierls ansatz on spin system histories, can be used effectively to compute local properties 
in stationary states when the underlying graph of interactions is sparse, and locally tree-like. A key technical step 
to make this a computationally feasible scheme is a stationarity assumption, here termed time factorization. With 
this step we find a message-passing scheme similar to Belief Propagation, but for dynamics, and for non-equilibrium 
systems. When these approximations work, at sufficiently high temperatures (weak interactions) they offer a huge 
speed-up compared to Monte Carlo 

We now give a synthetic overview of the dynamic cavity approach. A total history of a collection of spins, evolving in 
discrete time steps according to some dynamics, can be described by the total joint probability p(a{Q) . . . <j(t)). From 
such a total joint probability one can construct the marginal probability over the history of one spin Pi(cri(0) . . .ai{t)). 
If the underlying graph of interactions is locally tree-like, then this marginal probability can be expressed using 
"messages" from the neighboring nodes of the type ^j^i{aj(0) . . .aj{t)), and these messages in turn obey update 
rules similar (in principle) to Belief Propagation. In contrast to cavity method in equilibrium analysis, messages 
carry the whole information of spin histories over time. The difficulty arises when we want to marginalize further to 
the configuration of one spin at one time in its history. In general, the resulting equations are non-Markovian, and 
hard to solve. One level of approximation is to assume that the messages factorize in time, as iij^i{aj(0) . . . (yjit)) = 
Y[l=o ^'^d stationary state we can further assume that the terms in the product (^^!^j(crj(s)) are all 

the same function i.e. do not depend on s). The terms in the factorized messages {jij^iiaj)) then obey a set of 
distributed equations analogous to (but more complex than) Belief Propagation. Depending on the dynamics of the 
spin system, the resulting equations will differ. We here investigate both parallel updates, as studied recently by Neri 
and BoUe [22: as well as random sequential updates. 
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A summary of our results is as follows. We study stationary states in a diluted asymmetric (kinetic) Ising model. 
This well-studied non-equilibrium model has (qualitatively) three parameters: the asymmetry degree (how much the 
system is out of equilibrium) ; the connectivity of the underlying interaction graph (how accurate a Bethe-Peierls ansatz 
can be expected to be), and the strength of the interactions, customarily denoted inverse "temperature" . Furthermore, 
the type of the dynamics influences behavior. We here study parallel (synchronous) updates and sequential updates 
(one spin at a time). We compare dynamic cavity results to numerical simulations of the spin system dynamics, which 
we will refer to as Monte Carlo. The first result is that the more asymmetric the network, the better dynamic cavity 
method agrees with Monte Carlo. This can be understood as an effect of lack of memory in the graph: message go out, 
and rarely come back. The second result is that the sparser the underlying graph of interactions, the better dynamic 
cavity method agrees with Monte Carlo. This can be understood as an effect that the Bcthe-Peierls approximation in 
itself being more accurate for locally tree-like graphs. The third result is that dynamic cavity agrees well with Monte 
Carlo at high temperature, but deviates from numerical simulations of the full dynamics at low enough temperatures. 
The way in which dynamic cavity diverges from Monte Carlo at low temperature depends on the update rule. For 
random sequential dynamics we find that dynamic cavity (in the time factorized approximation for the stationary 
state) goes to a fixed point, but which does not agree with the stationary state estimated from Monte Carlo. For 
parallel updates on the other hand, we find that dynamic cavity method at low enough temperature does not go to 
a fixed point. As the system size increases, we find that dynamic cavity matches Monte Carlo better, and that the 
fluctuations in parallel updates diminish. In summary, the parameter ranges where dynamic cavity method agrees well 
with Monte Carlo seem to be rather wide, and not decreasing with system size. We therefore believe that dynamic 
cavity as presented here will be a promising avenue to compute and explore stationary states of large non-equilibrium 
systems. 

The paper is organized as follows. In section( [ll| we briefly describe the model, introduce the macroscopic observables 
and we study the dynamic update rules. Section! Ill I summarizes previous studies using "naive" mean field theory for 
the kinetic Ising model. In section(IV) we introduce the dynamic cavity method for parallel and sequential updates, 
and in section ([v|) we compare dynamic cavity to direct numerical simulation (Monte Carlo). Section! VI I contains 
concluding remarks. 



II. THE DILUTED ASYMMETRIC (KINETIC) SPIN GLASS 

Kinetic Ising models were originally motivated by neural networks, to extend the Hopfield model to asymmetric 
interactions |13j . These non-equilibrium systems with random interactions have formal similarities to the equilibrium 
and especially out-of-equilibrium (relaxation) dynamics of spin glasses, and therefore have a long history of study using 
methods from that field. Sompolinsky and Zippelius f31' introduced a formalism based on the Langvin equations of 
spherical spin models. An analogous approach was then proposed by Sommers in |30J through a path integral analysis 
of the Glauber dynamics. More recently, the dynamic replica theory has been developed, partly with an application 
to this kind of systems in mind 01 1201 121| . A general feature of dynamic replica theory is an average over disorder 
(average over a class of random graphs and random interactions) ; in addition technical assumptions may be needed 
such as considering the dynamics of some finite set of appropriate observables [3] . In [TD] a combination of dynamic 
replica theory and the cavity method (equilibrium) concept was applied to finitely connected disordered spin systems. 
An alternative approach to dynamic replica theory is generating functional analysis (GFA). The GFA has a long 
history in analyzing non-equilibrium statistical mechanics of disordered systems |8]. In particular, it allows us to 
study systems with non symmetric interacting couplings |TT] . More recently it has been developed to study the 
dynamics of spin glass systems with finite connectivity interactions In its general formalism, GFA aims at 

solving the dynamics of spin system exactly, however, due to the complicated nature of problem, one needs either 
to use a perturbative approach [16J or to restrict the analysis of GFA up to some approximation levels [TT1I21]. For 
ferromagnetic systems with regular connectivity, which is much simpler than spin glasses, a recursive set of dynamical 
equations can be derived for some finite macroscopic observations [29| . 

The dynamic cavity method shares some of the features of the dynamic replica theory and generating functional 
analysis. In equilibrium analysis, the Ising spin glass systems in Bethe lattice model have been solved by cavity 
method on the level of replica symmetry breaking. However in contrast to replica method, cavity applies to one 
single graph instance (one single set of interactions). Neri and Bolle [33] and Kanoria and Montanari [T3| considered 
dynamic cavity method in parallel updates under respectively Glauber dynamics and majority rule dynamics. In the 
present work we extend the dynamic cavity method to random sequential updates, and investigate stationary states of 
the diluted asymmetric spin glass over the wide range of parameters, for both parallel updates and random sequential 
updates. 

The asymmetric diluted Ising model is defined over a set of N binary variables a = {(Ti, . . . , cn}, and an asymmetric 
graph G — (V, E) where y is a set of N vertices, and E' is a set of directed edges. To each vertex Vi is associated 
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a binary variable cr,. The graphs G are taken from random graph ensembles with bounded average connectivity. 
Following the parameterization of we introduce a connectivity matrix Cy-, where Cij = 1 if there is a link from 
vertex i to vertex j, Cij — otherwise, and matrix elements Cij and Cki are independent unless {kV\ — {ji}. The 
random graph is specified by marginal (one-link) distributions 

c c 

and the conditional distributions 



p{c^J I Cji) = eSc^^^c,, + (1 - e)p(cij) . (2) 

In this model the average degree distribution is given by c, and the asymmetry is controlled by e G [0, 1]. The two 
extreme values of e give respectively an asymmetric network (e = 0) , where the probabilities of having two directed 
links between pairs of variables are uncorrelated, and the fully symmetric network (e = 1) where the two links i ^ j 
and j ^ i are present or absent together. The parameter set is completed by a (real-valued) interaction matrix 
We will always take Jij to be independent identically distributed random variables with zero mean and unit 
variance (Gaussian or binary) such that for the fully connected networks (c = N), the the interactions scale as the 
Sherrington-Kirkpatrick model. 

The definition of the kinetic Ising model is completed by a dynamics, or a spin update rule. In the synchronous 
update rules, which will be considered here, at each (discrete) time, a set of candidate spins are selected, and then 
updated according to the rule 

. +1 with probability {1 + exp(2/3 h.it + At))}-^ 

a,[i + i2.i) I _i .^vith probability {l-fexp(-2/3/i,(t + Ai))}-i ^"^^ 

where At is the time interval in which the update procedure takes place and hi{t) is the effective field on spin i at 
time step t 

h^it) = J2 ^-f<Jj{t~At) + e,{t) . (4) 

and the parameter /?, analogous to inverse temperature, is a measure of the overall strength of the interactions. The 
notation j e di in ([s]) and (|4| indicates all vertices having a direct links to node i (defined by Cji = 1) and 9i is the 
(possibly time-dependent) external field acting on spin i. 

We will consider two cases of synchronous updates: either all spins are selected and updated at each time step, or 
only one spin is randomly selected and updated in each time step. We refer to the first update rule as parallel, and 
the second as sequential. The time interval between updates is taken At = 1 in parallel updates, and At = in 
sequential updates, such that in both cases 0{N) spins are updated per unit time. 

The joint probability distribution over all the spin histories p((t(0), . . . , (j{t)) has in principle the following simple 
Markovian form 

t 

Pirn . . . , a{t)) ^ n ^[^(^) I h{s)]p{m) (5) 

where W is the appropriate transition matrix describing dynamics and updates. Solving equation ([s]) is in general 
infeasible for large system size, since it consists of 2* 2^ equations corresponding to all possible spin history configura- 
tions. One therefore needs to restrict the analysis to some restricted set of observables. In this paper we are interested 
in approximations built on marginal probabilities. The evolution of a a single spin is (trivially) defined by summing 
over the histories of all spins except one 

p,(a,(0),...,a,(t)) = P(f?(0),...,a(t)) (6) 

5'\,(0),...,ffy,(t) 

and similarly for pairwise joint probability of the histories of two spins (t > t') 



K,(a,(0),...a,(t),a,(0),...,fT,(t')) - p(f?(0), . . . , a(t)) 

CT\,,j(0),...,ff\i.j(t) 



(7) 
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Consequently, the time evolution of single site magnetization and the pairwise correlations can be obtained from EqQ 
and Eq0 

K(cri(ljj = ^ (Oj 

p,j(fT,(t),aj(t')) = [l + TO»Wa.(t)+TOj(t')CTj(*') + Czj(t,t')'^»Wf^j(i')]/4 (9) 

Computing the marginal probabilities directly is clearly an intractable problem for large system sizes, since the exact 
enumeration requires summation over an exponentially large number of states. The situation is even more involved 
for time dependent quantities since in addition the effect of past history must be taken into account. 



III. NAIVE MEAN-FIELD APPROXIMATION 



The dynamics of spin glass models has been widely studied using mean-field approximations; we follow recent 
practice in referring to this level of approximation as naive mean-field j25| . In 5, 6] such a theory was proposed to 
describe the dynamics of Little-Hopfield model in the case of fully asymmetric networks. Here we briefly summarize 
the naive mean field theory for diluted asymmetric systems. The time evolution of magnetization and correlations for 
the parallel update defined in Eq([5| can be explicitly written as 

m,(t) = (tanh(^/i,(t))) (10) 

where (...) represents the average over probability distribution at time t. Similarly for the sequential update we have 

m,(t+^)-(l-^)m,(i) + ^ (tanh(/3/i,(t))) (11) 



Equations (10 and (11) are yet exact. The right hand side of equations is however not easy to compute. The 
assumption in naive mean field theory is to substitute the effective field hi (t) at time t with a time dependent random 
Gaussian noise which does not contain spins configurations [5] . The formula for the time evolution of the magnetization 
(expectation value of a single spin) is then 

TOi(i+l) = tanh (^/3hi{t)j parallel update (12) 

mi{t + —) — mi{t) + — (^ta,nh{l3 hi{t)) ~ m,i{t)^ sequential update (13) 

where hi (t) is given by the mean values of spins neighboring i and a Gaussian noise reflecting the effect of neighbors on 
the dynamic of spin i (see [S]and |30|V The fixed point of two dynamic updates coincide, but system size is important 
for the stability of fixed point in the sequential update whereas it has no effect on the parallel update scheme [2 [32] . 
The naive mean-field approximation, although introduced quite some time ago, remains a main theoretical tool to 
analyze kinetic Ising models. More recently, these approximations have been used as the basis for "kinetic Ising" 
reconstruction schemes [HI [IHl [31] ■ 



IV. THE DYNAMIC CAVITY METHOD 



In this paper we use the terms Bethe-Peierls approximation (BP) and cavity method interchangeably. Their modern 
use grow out of replica theory in spin glasses, but in a form which can be applied to a single instance [H [T7| [T^] . We 
use the term dynamic cavity method for the use of the Bethe-Peierls approximation on spin histories. 

The main idea of BP is to ignore long loop correlations, since BP is exact on trees. BP can then be described by 
a set of auxiliary graphs called cavity graphs, which are identical to the original graphs but with one of the vertices 
and its associated variable removed. The effects of removing a target node i appear in the effective fields Eq([4]) acting 
on the set of variables neighboring node i with direct interactions outgoing from spin i. Fi^ illustrates the argument 
graphically. By the assumption that short loops are absent and long loops are ignored, the variables associated to the 
vertices which were connected to the removed vertex (denoted by i in Fi^l]) are independent in the cavity graph. We 
can then consider the joint probability of the histories of all the other spins p^'\(7\i{0) , ■ ■ ■ , (f\iit) |^'^''(0), . . . , Ol^^it)) 

under the influence of the external fields modified by the action of spin i i.e. 9j^\s) — Oj{s) + Jjiai{s — At). Since 
the removed vertex only modifies the effective fields of its neighbors with an outgoing edge we can further marginalize 



FIG. 1: Left figure: a directed network representation of the asymmetric Ising model. Interacting pairs are connected by 
directed edges where the arrows indicate the the asymmetric nature of model. Right figure: the cavity graph created by fixing 
the spin i to have value Si. This effects other spins which were connected to vertex i by an incoming edge from i (shown by 
red color). In graphs with tree structure, this leads immediately to a factorized probability distribution for the set of spins in 
different branches outgoing from node i. 



the joint probability distribution to the set di of spins "in the cavity" which were directly connected to vertex i, 
P^^H'^diiO), . • • , (Jdi{t) l^ai(0)j • • ■ I (*))■ The assumption of no loops means that the set of spins neighboring vertex 
i "in the cavity" produced by removing vertex i are independent, and that therefore this joint probability over these 
spin histories is factorized as 

p«(aa.(0), . . . |^^^^(0), . . . , e^{t)) = n M,^.(^,(0), . . . , o,(t)\ef(fi), . . . (t)) (14) 

The above approximation is exact in trees, but in general it can be only considered as an approximation. For random 
graph ensembles with diluted interactions, the typical loop length diverges in thermodynamic limit and BP is therefore 
expected to become accurate. 

By the same argument we can marginalize the joint probability distribution over any single vertex in the neighbor- 
hood of i and consequently all individual spins. The cavity assumption imposes the marginal probability fii^j to be 
dependent on the spins that are directly connected to the vertex j. Therefore we can interpret marginal probabilities 
in the cavity graph to be a set of "messages" exchanged among interacting pairs. These messages themselves obey 
recursion equations ( "BP update equations" ) which for parallel update are 



E 



n A*.^,K(o), ...,a,it - At)\ei'\o), ei'\t - At)) 



Ul=^wA'^M\hf{^'^))^^,^^{<Jm (15) 



Here /i^*^ is the effective field on spin j in the cavity graph 

Jk 



hf(s)= J2 ^<Jk{s-At) + 0,{t) 



(16) 



k£dj\i 



and WjicTj \ h^j^\s)) is the transition probability for the single spin j in the cavity graph. Similarly for the sequential 
update we have (see appendix for details) 



/i,-..(a,(O),...,a,(i)|0y^(O),...,0].^^(i))= E.,^^,io),....,.oMt-At)nke3J\^^^k^M^^ 

nl=i \7f^ji<^jis) I hf{s)) + (1 - ^)<5.,(«),.,(.-At)l H^^{'JJm (17) 
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Fig j2] illustrates how messages are distributed among interacting vertices: using the terminology of belief propagation 
(BP), the conditional probability fij^i{aj{0, . . . ,aj(t))) can be interpreted as a message sent from variable j to its 
neighbor i to indicate the probability of observing spin history crj(O) . . . <Jj{t) when i is removed from the network. 





FIG. 2: Dynamic message passing scheme for the diluted Ising systems. Each message contains information describing the 
evolution of marginal probability in the cavity graph. Messages are exchanged among the interacting pairs. The history of 
target vertex may effect the incoming message. 

The probability of the full history of one spin ( "BP output equation" ) is on the same level of approximation 



(Tai(0)...CTa,(t-At) 



At)) 



nl=i WMis)\h^{s))MMO)) 



(18) 



A peculiarity of this formula (shared with Eq( 15 ) and Eq( 17 )) is that the probability distribution of i depends on the 
neighbors di through the effective field hi{s), but the messages sent from the neighbors to i also depend parametrically 

(i) . . . 

on the history of i through the modified external fields 0f, . This difficulty is absent for fully asymmetric networks, 
since then the messages sent to i do not depend on the history of «, as we briefly review below in Subsection |IV A| In 
the general case of partially or fully symmetric networks this difficulty is addressed in long time limit evolution by a 
stationary assumption and the time factorization ansatz in Subsection |IV B 



St arti ng from a suitable initial conditions for messages, the evolution of messages can be followed from equation ( 15 ) 
and ( 17|. Each message contains 2* different states corresponding to the trajectory of vertices. Marginalizing equation 



(15) and (17) to the single time step t would provide us the time dependent magnetization through equation (|8|) 



J2 M,^.(^,(O),...,a,(t)|0f (0),...,f?f (i)) 

<Tj(0)...<Tj(t-l) 



(19) 



However due to the history of target vertex i the single time message does not obey a Markovian process. Therefore 
one can only hope to proceed with this iterative procedure for only a few initial time steps. 

A. Dynamic cavity method in fully asymmetric networks 

Fully asymmetric diluted Ising models where if spin i connects to spin j then spin j does not connect back to spin i 
have, as alluded to above, the simplifying property that influences (through interactions) do not return. Accordingly, 
the dynamic cavity models also simplify. In the model family considered here, e = corresponds to the case where the 
probabilities of having two directed links between spin pairs are independent. For fixed connectivity, the possibility of 
having simultaneously two links between spin pairs can be neglected in the thermodynamic limit, this case can hence 
be assimilated to a fully asymmetric network. For the dynamic cavity messages we then have 



^fc^i((Tfc(0), . . . , ak{t)\ai{0), . . . , crj(t)) = ^fe_>i((TA;(0), . . . , <7k{t)) 



(20) 
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In parallel update this property allow us to sum over all the history except the last time step from Eq( 15 ) and Eq( 19 1 
resulting in 



E n 

3oi\j{t-l) kedt\j 



1)) w.,{a,is)\h[=\s)) 



(21) 



where we have introduced ^l{ai{t)) as the message sent from i to j at time step t. Despite the general non-Markovian 
dynamics of the marginalization to one time instance in dynamic cavity, for parallel update the fully asymmetric 
case does follow a Markovian dynamics: at any time, the messages carry information from the incoming messages 
at only one time step before. Note that, resulting Markovian dynamics in this case is a consequence of update 
rule and asymmetry properties in the coupling interactions and is not necessarily restricted to belief propagation 
approximation. The exact dynamics in fully asymmetric case would follow a similar equation for the joint probability 
distribution in which at each time the information from one time step is required. However, computationally, this is 
unfeasible to study the evolution of large set of coupled spins. The case of random sequential update is more delicate, 
and will be discussed separately [T]. 

We now consider the transition matrix to represent the Glauber dynamics 



Wi{ai{s) I hi{s)) 



2 cosh(/3 /i,(s)) 

from which the single site magnetization under parallel updates follows 



E 

^a^{t-l) 



(22) 



(23) 



Due to the particular form of transition matrix in the Glauber dynamics (local normalization) , the resulting dynamic 
BP equations in the fully asymmetric networks yet requires a summation over the whole configuration of neighbor- 
ing vertices. Note that some more simplified transition matrices would provide us with an even more efficient BP 
approximation where the number of required summations scales linearly with the size of neighboring vertices TJ . 

As an interesting extension, the fully connected graphs with weak interactions can be realized by taking the limits 
c — iV and N ^ oo. Since the interaction couplings are scaled with 1/c the variables become weakly connected in this 
limit and the graph statistically uncorrelated. Introducing 5hi(t) = X^jeOi ^i'^ji^ ~ 1) ~ ^ji^ ~ 1)) can expand 
the right hand side around Shi considering the fluctuation of spins to be small with respect to their mean value at 
each time 



mi{t) = tanh 



/?(^^m,(t-l) + 0,:) 



0{c- 



(24) 



To first order in 1/c we end up with the same equation for naive mean field approximation as the one introduced 
by [3] and to the second order in 1/c with the dynamic TAP approximation [121 127] . 



B. Dynamic cavity method and stationary states 



In this subsection we introduce an approximation scheme for general (not fully asymmetric) networks assuming 
that the system is in a stationary state. The problem to be solved is the non-Markovian nature of the evolution 
of the single-spin single-time marginals {e.g. Eq(19)) which follows from marginalization of the full dynamic cavity 
equations (15) and (17) over time. The approximation is that the dynamic cavity messages factorize over time: 



M»^j(ct»(0), . . . ,CT»(i) I cTj(0), . . . , a jit)) = Yl f^Uji'^ds) kj(s - 1)) 



(25) 



This time-factorization assumption is clearly not appropriate to describe transients, where we would expect both 
dependencies between messages and that the functional form of the messages depend on time. However, in a stationary 
state it may be acceptable to take the messages independent in time, and it is reasonable to assume that the single- 
time messages do not depend explicitly on time. On the other hand, from the computational point of view, time- 
factorization provides us a closed set of equations for the single-time marginals, which makes the whole scheme 
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computationally feasible: 



In above, for parallel updates, as treated in [53], the summation over time history is resulted in n\l^j{cri{t — 2)). 
We refer to appendix for the corresponding equations for sequential updates. Note that in this approximation the 
single-time dynamic cavity messages at time t depend on messages sent at most two time steps earlier. For Ising spins, 
we can write the single-time dynamic cavity messages using cavity biases /^*_^.j(o'i(^)) = :fco'sh(i*"* ' ^^'^ inserting 



this equation into ( 26 1 we get an evolution equation for the cavity biases 



2/3 



E 

.■5aAj(*-l).'^'(*-2) 



Ylkai\j 2 cosh[/3 {uk^i + JikCTiit - 2))' 

g^«.^j(i-2)(T,(t-2) 



g/3h«>(t)a.(t) 



2cosh{l3h[^\t)) 2cosh{l3u,^j{t - 2)) 



(27) 



supplemented by Eq|l6| i.e. h^^\t) — J2k<£dj\i ~ 1) + ^ji^) the cavity fields. Solving for the stationary 

state of the kinetic Ismg mod el u sing dynamic cavity equations in the time-factorized approximation hence means to 
find a fixed point of (|27[) and ( 16 1 when the external fields 9 are independent in time. Note that ai{t — 2) contributes 



only when the edges from i to j and from j to i are both present. Therefore, in the fully asymmetric network this term 
disappears and we get back to Eq(21 1. On the other hand, for fully symmetric networks, as has been already pointed 
out in 



and also in GFA analysis in [3], the ordinary belief propagation equation is a solution of dynamic cavity 
equation in the time factorized approximation. Indeed, it can be verified that Eq( 27) for CijJij = CjiJji admits a 
solution of the form Ui^j ~ 0i + Pj^k^di^j tanh(/3Jfej tanh{uk^i)) that is the ordinary belief propagation equation 
for Ising systems with pairwise interactions. We note also that in a transient we can compute the time evolution of 
magnetization from Eq([6]) 



rrii (t) 



E 



n 



kdi\j 



2cosh[/3 {uk^i + Jtk(^i{t - 2))] 



tanh 



/3(^ Jj,aj(i-l) + 0, 

jedi 



2cosh(/3u,^j (t- 2)) 



(28) 



This is not expected to be accurate unless we are already in a stationary state, but will be used below in Section [V] 
to monitor the approach to the stationary state. 



V. RESULTS 



In this section we investigate the stationary states of diluted spin glass by dynamic cavity approach in both parallel 
and sequential update. We show how the total magnetization of diluted Ising systems as computed by dynamic cavity 
method evolves with time. In order to verify the results, we perform a numerical simulations (Monte Carlo) based on 
the appropriate dynamics. 



A. BP versus Glauber dynamics 



We first generate diluted graphs from random ensemble of size 10'^ and 10^ using asymmetry and connectivity 
parameters as in EqQ and Eq([2| . For each graph instance we iterate the dynamic cavity equations and simulate Monte 
Carlo analysis to test the accuracy of dynamic cavity method. The system is initialized to a random configuration 
with small external fields acting on each individual variable. In Monte Carlo simulations, we generate up to 10000 
samples to estimate the time evolution of magnetization at short time steps. We focus on the high-temperature regime 
to avoid the spin glass phase (for networks with symmetric interaction couplings, and presumably also for weakly 
asymmetric networks). Typical results for the time evolution of total magnetization m{t) — '^/N^^^-ymi{t) under 
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sequential updates are illustrated in FigjSj For short times, dynamic cavity method differs significantly from Monte 
Carlo, and also displays its own (unrelated) dynamics. For long times however dynamic cavity reaches a fixed point, 
and this stationary state agrees well with Monte Carlo (for the magnetizations). The time needed for Monte Carlo 



0.015 - 




0.025 - 




0.03 




FIG. 3: Time evolution of magnetization in a diluted networks under the sequential update. The interactions are random 
variables from normal distribution, connectivity parameter is fixed to c = 3 and simulations are performed for e = 0, 0.5, 1 and 
/3 — 0.7, 2, 3. A non-zero static external fields 6i — 0.01 acts on each individual spin. Left panel, fully asymmetric network, 
middle panel partially asymmetric with e = 0.5, right panel a fully symmetric network. Symbols show the dynamic cavity 
method and solid lines show the numerical simulations by Monte Carlo averaged over 10'* samples of the size A*' — 1000. 

to reach the stationary state grows with /3 as it can be seen from the FigjS] while for dynamic cavity in sequential 
update there is no such noticeable dependence except for the fully symmetric networks (right panel). 

For parallel dynamics, the situation is more involved. For low enough /3, again we have consistent results for 
dynamic cavity and Monte Carlo simulations (Fig. |4j left panel). Note that the dynamic cavity method for diluted 
fully asymmetric network in parallel update is exact, and (Fig. |4j left panel) lower curve gives a measure of the 
numerical fluctuations. In (Fig. |4j right panel) we verify that for symmetric interactions and in low enough /3, the 

solution for dynamic cavity is also a solution of ordinary BP. Introducing D{t) = ^ / N ^^^^{rni(t) — rn^^^^Y where 

m|^^' is the single magnetization obtained from output belief propagation (equilibrium), we expect to get a zero 
value for D{t) by evolving dynamic cavity during time. Fig. [4j right panel, shows the evolution of D{t) for /3 = 1 
and (3 — 2. For both cases, ordinary BP converges to a fixed point and dynamic cavity method predicts stationary 
solution to the time dependent magnetization. For /3 = 1 a, fast convergence for dynamic cavity solution is observed 
to the total magnetization predicted by ordinary BP. Increasing /3 (still in the phase where ordinary BP converges) 
would require longer relaxation time for constant results. 

In the large /? limit, however, the system may fall into limit cycles with no stable stationary state. This can be 
observed by computing total magnetization over long periods of time. In Fig[5]the time evolution of total magnetization 
in a fully symmetric networks computed by dynamic cavity is illustrated for /3 = 4 and (3 = 5. We observe that even 
in long time limit, the dynamics according to cavity method, does not approach a fixed solution but roughly oscillates 
between two states. Interestingly, for fully asymmetric networks we do not observe such cyclic behavior. 



B. Cyclic stationary states 

In order to study the cyclic behavior of stationary states in parallel update we introduce a time dependent quantity 
which measures the difference of total magnetizations in two successive time steps 

1 ^ 

4=1 

A zero value for S{t) indicates the existence of stationary state. On the other hand, non-zero 5{t), even after long times, 
means that either a stationary state does not exist, or is not reachable in finite time. Our results (from Fig. [5] and 
data not shown) are that cyclic behavior (of this type) , is observed for parallel updates at sufficiently low temperature 
if the network is not fully asymmetric, but not for sequential updates. The strongest cyclic behavior belongs to fully 
symmetric networks. At fixed connectivity the amplitude of the oscillations tend to decrease with system size (data 
not shown). 
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FIG. 4; Left panel : the time evolution of total magnetization at /3 = 1 for networks of size TV = 10^ with asymmetric and fully 
symmetric interactions. In all cases, a stationary solution according to dynamic BP exist which is consistent with numerical 
simulations. Right panel: the comparison between fixed point of ordinary BP and the time evolution of magnetization obtained 
by dynamic cavity method for fully symmetric networks. Inverse temperature /3 is chosen to be 1 and 2. In both cases, ordinary 
BP converges and the limit D{t) — > exists. For lower /? a faster convergence of dynamic cavity method is found. 




FIG. 5: Dynamic cavity results for diluted spin glass systems in parallel update at low temperature. The interactions are 
random variables from normal distribution, the system size is A'' = 1000 with an average connectivity c = 3. Left panel, the 
time evolution of magnetization for e = 0.5 and /3 = 4, 5 in absent of external fields. It shows an oscillatory behavior even at 
long limit time (for non- fully asymmetric and low temperature). Right panel, the evolution of 6 for system size oi N = 1000 
after t = 10"* steps and for e = 0.2, 0.5, 1. For fully asymmetric networks 3 = showing the existence of stable stationary state. 

VI. CONCLUSION AND DISCUSSION 

In this paper we have studied stationary states of diluted asymmetric (kinetic) Ising models by the dynamic cavity 
method, both for parallel updates (as in and [T3]) and for sequential updates (new). We find that for a many such 
systems with different asymmetry, sparseness and interaction strengths, total magnetization computed by dynamic 
cavity matches direct numerical simulation (Monte Carlo) - at many orders of magnitude less computational effort. 

Nothing does however come for free, and for all cases except a fully asymmetric network in parallel update, the 
dynamic cavity method presumably fails at low enough temperature. We observe (for the cases we have studied) 
different failure modes for different update rules, where for parallel updates the method simply no longer converges 
(no longer finds a stationary state), while for sequential updates the method converges to fixed point, but which does 
not correspond to the stationary state found by direct simulation. 

The full phase diagram of these models remains to be done; we have e.g. not verified that dynamic cavity fails at 
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low temperature for all models except fully asymmetric ones (although we expect that to be true). The comparison 
should also be done in a more detailed spin-by-spin manner, and extended at least to pair-wise correlations. We 
intend to return to these topics in a future contribution. 

The sequential update model is close to an asynchronous update model (continuous time) i. e. the master equation, 
as simulated e.g. by the Gillespie algorithm. But at least for finite N, the two models are not identical, and the 
analysis carried out here should be extended to the master equation (if this is possible). 

Finally, since kinetic Ising models have recently been used for inference and dynamic network reconstruction [25 1 126 ) . 
and since ordinary BP has been used for inference in equilibrium systems [TS] it would be interesting if the two strands 
could be combined. 
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A. Appendix 

The dynamic cavity approach to the parallel update was first derived by Neri and Bolle in |23j . Here we show 
how it can be extended to sequential update. In the sequential update we assume that at each time step, one single 
variable is selected randomly and will be updated according to Eq( [s]) . The time evolution of probability distribution 
of all variables, follows 



Pirn, • • • , '^it)) - n ^('?(*) I ks))pim) 



(30) 



where the transition probability contains the sequential update. The appropriate choice for the single variable update 
is the uniform probability distribution over all spins (see [5]) 



(t-At) Wi{(7t(t) I hi{t)) 



(31) 



Following the cavity graph argument discussed in section IV the evolution of marginal probability follows 
Pi 



,(a,(o),...,a,(i)) = ^ J2 n i^k^r{'^kio),...,<Jkit~At)\el:\o),...,ei'\t-At)) 

t N ( I 

n { ]l^'^kis),a,{s~At)Mc']is)\'T'jis)) } P^(C^»(0)) 

s=Atj=l [k^] I 



(32) 



The summation over variable update inside the formula will contribute in two terms: either the index j is equal to the 
cavity variable i or is different. We can now introduce messages among neighboring variables carrying the marginal 
probability in the cavity graph. They fulfill the same recursive equation as parallel update 



M,^,(a,(O),...,a,(t)|0p'^(O), 



.,^F^(t)) 



J2 n M.^«(^'c(o), ■ ■ • , '^kit - st)\el:\o), e\:\t - At)) 

^ai\]iO),----.^ai\]it) kedi\j 
t 



n 

s=At 



^ w^{(ji{s) I hi{s)) + (1 - ^)(5<Ti(s),<Ti(s-At) 



N 



li,^ji<7,{0)) 



(33) 



which is Eq(17l used in the section(IVl for dynamic cavity in the sequential update. In the stationary state where the 



initial condition is assumed to be irrelevant we can perform the time-factorized approximation over the time evolution 
of messages. The time evolution of messages at each time then simplifies following a Markovian dynamics of length 
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two. 

^^UMit)\o?\t)) = - At)) + (i - ^) E IT ^ I - 

(Ti(t-2At),CTai\j(t-At) kedi\j 

w^{a,{t) I /,„(t)) M-^f - 2At)034) 
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